5261project

2. Descriptive Statistics

Means SDs Skewness Kurtosis Betas

	Mean (price)	Mean(return %)	SD (price)	SD(return %)	Skewness (price)	Skewness (return %)	Kurtosis (price)	Kurtosis (return %)	Beta
MSFT	117.508	2.495	86.817	5.918	1.043	0.231	-0.080	0.691	1.004
TSLA	203.202	4.802	298.245	17.296	1.823	1.293	1.974	2.342	1.993
AAPL	58.870	2.774	45.183	7.744	1.271	-0.233	0.292	-0.215	1.247
TWTR	34.921	0.506	14.229	14.535	0.650	0.416	0.057	0.446	0.869
AMZN	1550.228	2.595	1062.254	8.163	0.549	0.402	-1.051	0.727	1.175
FB	171.007	1.636	80.593	8.123	0.713	-0.335	-0.240	2.892	1.219
NFLX	264.865	2.569	176.639	11.741	0.454	0.472	-0.997	1.115	0.941
AAL	33.150	0.112	11.580	11.493	-0.356	-0.096	-0.899	0.284	1.559
DAL	42.931	0.840	8.403	9.510	-0.057	-0.187	-0.600	2.746	1.183
BAC	23.756	1.398	9.369	8.191	0.696	-0.159	-0.265	1.109	1.438
NVDA	64.032	5.146	74.026	11.778	1.699	0.071	2.263	0.489	1.452
WBD	29.026	0.106	6.327	11.332	1.093	0.854	1.452	1.192	1.244
S.P500	2770.638	1.036	799.620	4.002	0.966	-0.393	-0.052	1.433	1.000

Plots

Equity curve

Stationary Test

MSFT	0.1
TSLA	0.1
AAPL	0.1
TWTR	0.1
AMZN	0.1
FB	0.1
NFLX	0.1
AAL	0.1
DAL	0.1
BAC	0.1
NVDA	0.1
WBD	0.1
S&P500	0.1

Hist, Boxplot, qqplot

Distributions

Asset	t-distribution	Normal Distribution	GED
MSFT	-139.66549	-271.25573	-276.23424
TSLA	-38.58013	-80.34437	-69.20914
AAPL	-112.15671	-220.30961	-218.90939
TWTR	-50.63192	-97.29781	-94.90164
AMZN	-107.95913	-209.86197	-209.57404
FB	-111.98310	-209.20876	-217.11195
NFLX	-73.72148	-138.40822	-140.89449
AAL	-73.60630	-140.92136	-141.25665
DAL	-96.76172	-178.09420	-184.85237
BAC	-108.18352	-207.92323	-209.41384
NVDA	-71.69179	-136.35263	-137.27907
WBD	-77.72453	-151.06318	-153.27377

Sharpe’s Slope

Assets	Sharpe’s Slope
MSFT	0.3012400
TSLA	0.2364437
AAPL	0.2662491
TWTR	-0.0141628
AMZN	0.2307159
FB	0.1137256
NFLX	0.1581712
AAL	-0.0522252
DAL	0.0134156
BAC	0.0838003
NVDA	0.3764550
WBD	-0.0534724

M to Y

	Monthly Mean	Monthly SD	Annual Mean	Annual SD
MSFT	2.4946532	5.917747	29.935838	20.49968
TSLA	4.8015901	17.296293	57.619081	59.91612
AAPL	2.7737763	7.743821	33.285316	26.82538
TWTR	0.5061345	14.535005	6.073614	50.35073
AMZN	2.5953512	8.163112	31.144214	28.27785
FB	1.6357391	8.122604	19.628869	28.13753
NFLX	2.5691372	11.741366	30.829647	40.67329
AAL	0.1117527	11.493264	1.341033	39.81384
DAL	0.8395770	9.510284	10.074924	32.94459
BAC	1.3984216	8.191269	16.781059	28.37539
NVDA	5.1459413	11.778168	61.751296	40.80077
WBD	0.1060357	11.332111	1.272428	39.25558
S&P500	1.0358383	4.002103	12.430059	13.86369

Pairewise

Covariance Matrix

	MSFT	TSLA	AAPL	TWTR	AMZN	FB	NFLX	AAL	DAL	BAC	NVDA	WBD	S&P500
MSFT	0.0035	0.0038	0.0024	0.0005	0.0027	0.0019	0.0026	0.0020	0.0014	0.0021	0.0034	0.0020	0.0016
TSLA	0.0038	0.0299	0.0063	0.0036	0.0045	0.0041	0.0055	0.0032	0.0020	0.0034	0.0059	0.0043	0.0032
AAPL	0.0024	0.0063	0.0060	0.0024	0.0028	0.0030	0.0025	0.0024	0.0019	0.0016	0.0047	0.0019	0.0020
TWTR	0.0005	0.0036	0.0024	0.0211	0.0023	0.0029	0.0018	0.0012	0.0023	0.0024	0.0019	0.0050	0.0014
AMZN	0.0027	0.0045	0.0028	0.0023	0.0067	0.0032	0.0056	0.0015	0.0008	0.0020	0.0042	0.0017	0.0019
FB	0.0019	0.0041	0.0030	0.0029	0.0032	0.0066	0.0036	0.0019	0.0008	0.0023	0.0028	0.0021	0.0019
NFLX	0.0026	0.0055	0.0025	0.0018	0.0056	0.0036	0.0138	0.0013	0.0001	0.0018	0.0051	0.0021	0.0015
AAL	0.0020	0.0032	0.0024	0.0012	0.0015	0.0019	0.0013	0.0132	0.0086	0.0050	0.0032	0.0039	0.0025
DAL	0.0014	0.0020	0.0019	0.0023	0.0008	0.0008	0.0001	0.0086	0.0090	0.0039	0.0024	0.0036	0.0019
BAC	0.0021	0.0034	0.0016	0.0024	0.0020	0.0023	0.0018	0.0050	0.0039	0.0067	0.0027	0.0043	0.0023
NVDA	0.0034	0.0059	0.0047	0.0019	0.0042	0.0028	0.0051	0.0032	0.0024	0.0027	0.0139	0.0008	0.0023
WBD	0.0020	0.0043	0.0019	0.0050	0.0017	0.0021	0.0021	0.0039	0.0036	0.0043	0.0008	0.0128	0.0020
S&P500	0.0016	0.0032	0.0020	0.0014	0.0019	0.0019	0.0015	0.0025	0.0019	0.0023	0.0023	0.0020	0.0016

3. Portfolio Theory

With Short Sale

library(quadprog)
R = 100*return[,2:13]
mean_p <- apply(R,2,mean)
cov_p <- cov(R)
sd_vect_p <- sqrt(diag(cov_p))
# min(mean_p)
# max(mean_p)
### With shortsale
M_p = length(mean_p)
Amat_p <- cbind(rep(1,M_p),mean_p)
mu_P = seq(0.07, 5.4, length = 300)
# Target portfolio means for the expect portfolio return
sd_P = mu_P # set up storage for std dev's of portfolio returns
weights_p = matrix(0, nrow = 300, ncol = M_p) # storage for return
for (i in 1:length(mu_P)) { # find the optimal portfolios
  bvec_p <- c(1, mu_P[i])
  result_p = solve.QP(Dmat = 2 * cov_p, dvec = rep(0, M_p), Amat = Amat_p, 
                    bvec = bvec_p, meq = 2)
  sd_P[i] = sqrt(result_p$value)
  weights_p[i, ] = result_p$solution
}
plot(sd_P, mu_P, type = "l", xlim = c(0,15), ylim = c(0, 6), lty = 3, lwd = 2, main = "With Short Sale") 
# plot efficient frontier (and inefficient portfolios below the min var portfolio)
mufree_p = mean(return$`Treasury Bill 3 month (rf)`)# input value of risk-free interest rate
points(0, mufree_p, cex = 4, pch = "*") # show risk-free asset
sharpe_p = (mu_P - mufree_p) / sd_P # compute Sharpes ratios
ind_p = (sharpe_p == max(sharpe_p)) # Find maximum Sharpes ratio
#weights_p[ind_p,] # print the weights of the tangency portfolio
lines(c(0, 15), mufree_p + c(0, 15) * (mu_P[ind_p] - mufree_p) / sd_P[ind_p], lwd = 4, 
      lty = 1, col = "blue") # show line of optimal portfolios
points(sd_P[ind_p], mu_P[ind_p], cex = 4, pch = "*") # tangency portfolio
ind2_p = (sd_P == min(sd_P)) # find the minimum variance portfolio
points(sd_P[ind2_p], mu_P[ind2_p], cex = 2, pch = "+") # min var portfolio
ind3_p = (mu_P > mu_P[ind2_p])
lines(sd_P[ind3_p], mu_P[ind3_p], type = "l", xlim = c(0, 25), ylim = c(0,30),
      lwd = 3, col = 'red') # plot the efficient frontier
for(i in 1:12){
  text(sd_vect_p[i], mean_p[i],colnames(return[,i+1]), cex=0.8)
}

### MVP
(mvp_meanreturn <- mu_P[ind2_p])

## [1] 1.816957

(mvp_sd <- sd_P[ind2_p])

## [1] 4.980327

weights_mvp <- weights_p[ind2_p,]
weights_mvp <- t(data.frame(weights_mvp))
colnames(weights_mvp) <- colnames(return[2:13])
knitr::kable(round(weights_mvp*100,2))

	MSFT	TSLA	AAPL	TWTR	AMZN	FB	NFLX	AAL	DAL	BAC	NVDA	WBD
weights_mvp	50.17	-4.98	13.48	4.89	3.85	13.16	4.31	-9.16	21.82	8.07	-5.96	0.34

(mvp_meanreturn_ann <- mvp_meanreturn*12)

## [1] 21.80348

(mvp_sd_ann <- mvp_sd*sqrt(12))

## [1] 17.25236

### Efficient Portfolio Frontier
EPF_mean <- mu_P[ind3_p]
EPF_sd <- sd_P[ind3_p]

### Tangency Portfolio
(tan_meanreturn <- mu_P[ind_p])

## [1] 5.4

(tan_sd <- sd_P[ind_p])

## [1] 9.86121

(tan_var <- tan_sd^2)

## [1] 97.24347

(tan_sharpes <- (tan_meanreturn-mufree_p)/tan_sd)

## [1] 0.4753989

VaR&ES

### Tail dependence can be seen among the assets, therefore we can fit
### our portfolio with multivariate t-distribution.

## MVP VaR&ES
library(MASS)
alpha = 0.05
return_mvp <- rowSums(data.frame(
  weights_mvp[1] * return[, 2],
  weights_mvp[2] * return[, 3],
  weights_mvp[3] * return[, 4],
  weights_mvp[4] * return[, 5],
  weights_mvp[5] * return[, 6],
  weights_mvp[6] * return[, 7],
  weights_mvp[7] * return[, 8],
  weights_mvp[8] * return[, 9],
  weights_mvp[9] * return[, 10],
  weights_mvp[10] * return[, 11],
  weights_mvp[11] * return[, 12],
  weights_mvp[12] * return[, 13]))
fitt_mvp = fitdistr(return_mvp,"t")
param_mvp = as.numeric(fitt_mvp$estimate)
mean_mvpfit = param_mvp[1]
df_mvpfit = param_mvp[3]
sd_mvpfit = param_mvp[2] * sqrt((df_mvpfit) / (df_mvpfit - 2))
lambda_mvpfit = param_mvp[2]
qalpha_mvp = qt(alpha, df = df_mvpfit)
VaR_par_mvp = -100000 * (mean_mvpfit + lambda_mvpfit * qalpha_mvp)
es1_mvp = dt(qalpha_mvp, df = df_mvpfit) / (alpha)
es2_mvp=(df_mvpfit+qalpha_mvp^2)/(df_mvpfit-1)
es3_mvp=-mean_mvpfit+lambda_mvpfit*es1_mvp*es2_mvp
ES_par_mvp = 100000*es3_mvp
VaR_par_mvp

## [1] 6285.691

ES_par_mvp

## [1] 8957.941

## Asset VaR
S0 = 100000
qnalpha = qnorm(0.05)

### MSFT
q_msft = as.numeric(quantile(return$MSFT, alpha))
VAR_msft = -S0 * q_msft
#VAR_msft

### TSLA
fit_tsla <- fitdistr(return$TSLA, "normal")
param_tsla = as.numeric(fit_tsla$estimate)
mean_tsla = param_tsla[1]
sd_tsla = param_tsla[2]
VAR_tsla = -S0*(mean_tsla+qnalpha*sd_tsla)
#VAR_tsla

### AAPL
fit_aapl <- fitdistr(return$AAPL, "normal")
param_aapl = as.numeric(fit_aapl$estimate)
mean_aapl = param_aapl[1]
sd_aapl = param_aapl[2]
VAR_aapl = -S0*(mean_aapl+qnalpha*sd_aapl)
#VAR_aapl

### TWTR
fit_twtr <- fitdistr(return$TWTR, "normal")
param_twtr = as.numeric(fit_twtr$estimate)
mean_twtr = param_twtr[1]
sd_twtr = param_twtr[2]
VAR_twtr = -S0*(mean_twtr+qnalpha*sd_twtr)
#VAR_twtr

### AMZN
fit_amzn <- fitdistr(return$AMZN, "normal")
param_amzn = as.numeric(fit_amzn$estimate)
mean_amzn = param_amzn[1]
sd_amzn = param_amzn[2]
VAR_amzn = -S0*(mean_amzn+qnalpha*sd_amzn)
#VAR_amzn

### FB
q_fb = as.numeric(quantile(return$FB, alpha))
VAR_fb = -S0 * q_fb
#VAR_fb

### NFLX
q_nflx = as.numeric(quantile(return$NFLX, alpha))
VAR_nflx = -S0 * q_nflx
#VAR_nflx

### AAL
q_aal = as.numeric(quantile(return$AAL, alpha))
VAR_aal = -S0 * q_aal
#VAR_aal

### DAL
q_dal = as.numeric(quantile(return$DAL, alpha))
VAR_dal = -S0 * q_dal
#VAR_dal

### BAC
q_bac = as.numeric(quantile(return$BAC, alpha))
VAR_bac = -S0 * q_bac
#VAR_bac

### NVDA
q_nvda = as.numeric(quantile(return$NVDA, alpha))
VAR_nvda = -S0 * q_nvda
#VAR_nvda

### WBD
q_wbd = as.numeric(quantile(return$WBD, alpha))
VAR_wbd = -S0 * q_wbd
#VAR_wbd

VAR_asset <- c(VAR_msft, VAR_tsla, VAR_aapl, VAR_twtr, VAR_amzn, VAR_fb, VAR_nflx,
               VAR_aal, VAR_dal, VAR_bac, VAR_nvda, VAR_wbd)
var_df <- cbind(names, VAR_asset)
var_df <- rbind(var_df,VaR_par_mvp)
var_df[13,1] <- c("MVP")
var_df[13,2] <- 6285.69101812641
colnames(var_df) <- c("Assets", "Var")

knitr::kable(var_df)

	Assets	Var
	MSFT	6918.17869693084
	TSLA	23502.7562034031
	AAPL	9898.52221650451
	TWTR	23279.5289061133
	AMZN	10763.0915501661
	FB	10561.0446590248
	NFLX	13072.710950181
	AAL	15706.0583956349
	DAL	12891.9240002285
	BAC	13041.1530377474
	NVDA	11974.7969532429
	WBD	14903.9099722124
VaR_par_mvp	MVP	6285.69101812641

Assets’ Sharpe’s Ratios

names_sh <- data.frame(colnames(return[1,2:14]))
sharpes_list <- rep(NA, 13)
for(i in 2:14){
  sharpes_list[i-1] = (mean(unlist(return[,i]))-mean(unlist(return[,15]))/100)/sd(unlist(return[,i]))
}
sharpes_list <- data.frame(sharpes_list)
shar_df <- cbind(names_sh,sharpes_list)
shar_df

##    colnames.return.1..2.14.. sharpes_list
## 1                       MSFT   0.30124000
## 2                       TSLA   0.23644366
## 3                       AAPL   0.26624907
## 4                       TWTR  -0.01416282
## 5                       AMZN   0.23071594
## 6                         FB   0.11372558
## 7                       NFLX   0.15817121
## 8                        AAL  -0.05222523
## 9                        DAL   0.01341556
## 10                       BAC   0.08380025
## 11                      NVDA   0.37645499
## 12                       WBD  -0.05347243
## 13                    S&P500   0.08091924

Without shortshale

### Without shortsale
#R = 100*return[,2:13]
#mean_p <- apply(R,2,mean)
#cov_p <- cov(R)
#sd_vect_p <- sqrt(diag(cov_p))
### With shortsale
#M_p = length(mean_p)
Amat_p_noss <- cbind(rep(1,M_p),mean_p, diag(1,nrow=M_p))
mu_P_noss = seq(min(mean_p)+0.0001, max(mean_p)-0.0001, length = 300)
# Target portfolio means for the expect portfolio return
sd_P_noss = mu_P_noss # set up storage for std dev's of portfolio returns
weights_p_noss = matrix(0, nrow = 300, ncol = M_p) # storage for return
for (i in 1:length(mu_P_noss)) { # find the optimal portfolios
  bvec_p_noss <- c(1, mu_P_noss[i], rep(0,M_p))
  result_noss = solve.QP(Dmat = 2 * cov_p, dvec = rep(0, M_p), Amat = Amat_p_noss, 
                    bvec = bvec_p_noss, meq = 2)
  sd_P_noss[i] = sqrt(result_noss$value)
  weights_p_noss[i, ] = result_noss$solution
}
plot(sd_P_noss, mu_P_noss, type = "l", lty = 3, 
     lwd = 2, xlim = c(0,15), ylim = c(0,7), main = "Without Short Sale") 
# plot efficient frontier (and inefficient portfolios below the min var portfolio)
#mufree_p = mean(return$`Treasury Bill 3 month (rf)`) # input value of risk-free interest rate
points(0, mufree_p, cex = 4, pch = "*") # show risk-free asset
sharpe_p_noss = (mu_P_noss - mufree_p) / sd_P_noss # compute Sharpes ratios
ind_p_noss = (sharpe_p_noss == max(sharpe_p_noss)) # Find maximum Sharpes ratio
#weights_p[ind_p,] # print the weights of the tangency portfolio
lines(c(0, 15), mufree_p + c(0, 15) * (mu_P_noss[ind_p_noss] - mufree_p) / sd_P_noss[ind_p_noss], lwd = 4, 
      lty = 1, col = "blue") # show line of optimal portfolios
points(sd_P_noss[ind_p_noss], mu_P_noss[ind_p_noss], cex = 4, pch = "*") # tangency portfolio
ind2_p_noss = (sd_P_noss == min(sd_P_noss)) # find the minimum variance portfolio
points(sd_P_noss[ind2_p_noss], mu_P_noss[ind2_p_noss], cex = 2, pch = "+") # min var portfolio
ind3_p_noss = (mu_P_noss > mu_P_noss[ind2_p_noss])
lines(sd_P_noss[ind3_p_noss], mu_P_noss[ind3_p_noss], type = "l",
      lwd = 3, col = 'red') # plot the efficient frontier
for(i in 1:12){
  text(sd_vect_p[i], mean_p[i],colnames(return[,i+1]), cex=0.8)
}

### MVP
(mvp_meanreturn_noss <- mu_P_noss[ind2_p_noss])

## [1] 1.977063

(mvp_sd_noss <- sd_P_noss[ind2_p_noss])

## [1] 5.110456

weights_mvp_noss <- weights_p_noss[ind2_p_noss,]
weights_mvp_noss <- t(data.frame(weights_mvp_noss))
colnames(weights_mvp_noss) <- colnames(return[2:13])
weights_mvp_noss

##                      MSFT          TSLA       AAPL       TWTR       AMZN
## weights_mvp_noss 0.487966 -1.764314e-18 0.06592605 0.05609903 0.03247623
##                         FB       NFLX           AAL      DAL        BAC
## weights_mvp_noss 0.1346882 0.02256666 -2.931383e-17 0.146045 0.04453321
##                          NVDA         WBD
## weights_mvp_noss 6.441445e-18 0.009699598

(mvp_meanreturn_ann_noss <- mvp_meanreturn_noss*12)

## [1] 23.72476

(mvp_sd_ann_noss <- mvp_sd_noss*sqrt(12))

## [1] 17.70314

sum(weights_mvp_noss)

## [1] 1

rownames(weights_mvp_noss) <- c("weight")
knitr::kable(round(weights_mvp_noss*100,4))

	MSFT	TSLA	AAPL	TWTR	AMZN	FB	NFLX	AAL	DAL	BAC	NVDA	WBD
weight	48.7966	0	6.5926	5.6099	3.2476	13.4688	2.2567	0	14.6045	4.4533	0	0.97

### Efficient Portfolio Frontier
EPF_mean_noss <- mu_P_noss[ind3_p_noss]
EPF_sd_noss <- sd_P_noss[ind3_p_noss]

### Tangency Portfolio
(tan_meanreturn_noss <- mu_P_noss[ind_p_noss])

## [1] 3.999688

(tan_sd_noss <- sd_P_noss[ind_p_noss])

## [1] 7.96967

(tan_var_noss <- tan_sd_noss^2)

## [1] 63.51564

(tan_sharpes_noss <- (tan_meanreturn_noss-mufree_p)/tan_sd_noss)

## [1] 0.412526

alpha = 0.05
return_mvp_noss <- rowSums(data.frame(
  weights_mvp_noss[1] * return[, 2],
  weights_mvp_noss[2] * return[, 3],
  weights_mvp_noss[3] * return[, 4],
  weights_mvp_noss[4] * return[, 5],
  weights_mvp_noss[5] * return[, 6],
  weights_mvp_noss[6] * return[, 7],
  weights_mvp_noss[7] * return[, 8],
  weights_mvp_noss[8] * return[, 9],
  weights_mvp_noss[9] * return[, 10],
  weights_mvp_noss[10] * return[, 11],
  weights_mvp_noss[11] * return[, 12],
  weights_mvp_noss[12] * return[, 13]))
fitt_mvp_noss = fitdistr(return_mvp_noss,"t")
param_mvp_noss = as.numeric(fitt_mvp_noss$estimate)
mean_mvpfit_noss = param_mvp_noss[1]
df_mvpfit_noss = param_mvp_noss[3]
sd_mvpfit_noss = param_mvp_noss[2] * sqrt((df_mvpfit_noss) / (df_mvpfit_noss - 2))
lambda_mvpfit_noss = param_mvp_noss[2]
qalpha_mvp_noss = qt(alpha, df = df_mvpfit_noss)
VaR_par_mvp_noss = -100000 * (mean_mvpfit_noss + lambda_mvpfit_noss * qalpha_mvp_noss)
es1_mvp_noss = dt(qalpha_mvp_noss, df = df_mvpfit_noss) / (alpha)
es2_mvp_noss=(df_mvpfit_noss+qalpha_mvp_noss^2)/(df_mvpfit_noss-1)
es3_mvp_noss=-mean_mvpfit_noss+lambda_mvpfit_noss*es1_mvp_noss*es2_mvp_noss
ES_par_mvp_noss = 100000*es3_mvp_noss
VaR_par_mvp_noss

## [1] 6291.464

ES_par_mvp_noss

## [1] 9034.092

##4. Asset Allocation

Efficient Portfolio

#efficient portfolio
new_dat = return[,2:13]
return_year = 12*colMeans(new_dat)

cov_mat = cov(new_dat)
eff_port = efficient.portfolio(er = return_year,cov_mat,target.return = 0.06,shorts = FALSE)
eff_port$er/12

## [1] 0.004999994

#check return
ret=0
for (i in 1:12)
{
  val = as.numeric(return_year[i])*as.numeric(eff_port$weights[i])
  ret = ret+val
  
}
ret

## [1] 0.05999993

#monthly risk
eff_risk = eff_port$sd/sqrt(12)
eff_risk

## [1] 0.02134934

# vaR = -S0*(u + qnorm(0.05)*sd)
S0 = 100000
q = qnorm(0.05)
eff_var = -S0*(eff_port$er/12 + q*eff_risk)
eff_var

## [1] 3011.654

#excel
cut_off = qnorm(0.05,mean = (1+eff_port$er/12)*S0,sd = S0*eff_risk)
var2 = S0-cut_off
var2

## [1] 3011.654

var2/S0

## [1] 0.03011654

#Expected Shortfall
S0*((-eff_port$er/12) + (dnorm(qnorm(0.05))/0.05)* eff_risk)

## [1] 3903.755

Tangency portfolio

risk_free = colMeans(return[,15]/100)

tan_port = tangency.portfolio(er=return_year, cov.mat=cov_mat, risk.free=risk_free,shorts=FALSE)
tan_port

## Call:
## tangency.portfolio(er = return_year, cov.mat = cov_mat, risk.free = risk_free, 
##     shorts = FALSE)
## 
## Portfolio expected return:     0.4214667 
## Portfolio standard deviation:  0.06885209 
## Portfolio weights:
##   MSFT   TSLA   AAPL   TWTR   AMZN     FB   NFLX    AAL    DAL    BAC   NVDA 
## 0.4830 0.0686 0.0903 0.0000 0.0452 0.0000 0.0000 0.0000 0.0000 0.0000 0.3129 
##    WBD 
## 0.0000

# w1 for tangency portfolio, w1 for risk free
# w1+w2 = 1
# 0.4799845 * w1 + 0.08543894 * w2 = 0.06

# 1.3041868% for the tangency portfolio and -0.3041868% for the risk free

A = matrix(c(1, tan_port$er, 1, as.numeric(risk_free)), ncol=2)
B = matrix(c(1, 0.06),ncol = 1)
solve(A)%*%B

##           [,1]
## [1,] 0.1276228
## [2,] 0.8723772

#vaR
S0 = 100000
q = qnorm(0.05)
eff_var = -S0*(0.05/12 + q*0.06885209/sqrt(12))
eff_var

## [1] 2852.626

eff_var/S0

## [1] 0.02852626

#ES 
S0*((-0.05/12) + (dnorm(qnorm(0.05))/0.05)* 0.06885209/sqrt(12))

## [1] 3683.158

5. PCA

correlation matrix of the returns on your assets.

AAL and DAL are the assets that most highly correlated with correlation = 0.78. If we only consider about the 12 assets, then NFLX and DAL has the lowest correlation = 0.0089. The NFLX also has a correlation with Treasury Bill 3 month (risk free) = 0.006.

Diversification will reduce risk for the assets that are negatively correlated. For example, if the correlation between 2 assets is -1, which is the perfect value for diversification, if one asset goes up 5%, the other asset will goes down 5%. So combining these 2 assets would minimize our risk to 0. In our 12 assets, we did not see any negative correlation. As a result, diversification will not reduce risk.

(note: here are a few assets that are slightly correlated such as NFLX and DAL, NFLX and AAL, TWTR and AAL, NVDA and WBD, MSFT and TWTR.)

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:fBasics':
## 
##     densityPlot

df <- return[,2:13]
# correlation matrix 
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) {
    usr <- par("usr")
    on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    Cor <- abs(cor(x, y)) 
    txt <- paste0(prefix, format(c(Cor, 0.123456789), digits = digits)[1])
    if(missing(cex.cor)) {
        cex.cor <- 0.4 / strwidth(txt)
    }
    text(0.5, 0.5, txt,
         cex = 1 + cex.cor * Cor) 
}

pairs(df,
      upper.panel = panel.cor,    # Correlation panel
      lower.panel = panel.smooth) # Smoothed regression lines

Parallel coordinate plot

parcoords(df, rownames = F, reorderable = T , queue = T, 
          withD3 = TRUE, alpha=0.5)

PCA Analysis

Another way to dig into the correlation between the return of the assets, we looked at PCA and its plots. PCA helps to reduce the dimension of the dataset and compute principal components that contain most of the information in the data. From the PCA Analysis, we found that first eight components captures almost 90% of the total variance, while nine components could explain 93% of the variance. The Scree Plot shows that the unexplained part drop quickly by adding PC1 to PC6, and gradually slower after that. Therefore, using around 6 to 9 PCs would explain the majority of the variance in our data.

To further developed the correlations, we also draw a biplot. According to the principle of the biplots, the variables are positively correlated if their vectors form an angle smaller than 90 degrees. From our plot, all the assets are having a acute angle, so we conclude that they all having a positively relationship.

df = return[,2:13]
model_pca = prcomp(df, center = TRUE, scale. = TRUE)
su = summary(model_pca)
su

## Importance of components:
##                           PC1    PC2     PC3     PC4    PC5     PC6     PC7
## Standard deviation     2.0761 1.3586 1.06307 0.93727 0.8879 0.80884 0.78736
## Proportion of Variance 0.3592 0.1538 0.09418 0.07321 0.0657 0.05452 0.05166
## Cumulative Proportion  0.3592 0.5130 0.60717 0.68038 0.7461 0.80060 0.85227
##                            PC8     PC9    PC10    PC11    PC12
## Standard deviation     0.70899 0.67221 0.56977 0.54539 0.44294
## Proportion of Variance 0.04189 0.03766 0.02705 0.02479 0.01635
## Cumulative Proportion  0.89415 0.93181 0.95886 0.98365 1.00000

screeplot(model_pca, type = "line", main = "Scree plot")

autoplot(model_pca,  loadings = TRUE, loadings.colour = 'blue',
         loadings.label = TRUE)

#model_pca$rotation[,1:2]

factor analysis

X = scale(df,center=T,scale=T)
n=fa.parallel(X)

## Parallel analysis suggests that the number of factors =  2  and the number of components =  2

df.fa = fa(X,nfactors = 4,fm="mle",rotate = "varimax")
fa.diagram(df.fa)

faa1 = factanal(df,factor=3,rotation = "varimax")
faa1

## 
## Call:
## factanal(x = df, factors = 3, rotation = "varimax")
## 
## Uniquenesses:
##  MSFT  TSLA  AAPL  TWTR  AMZN    FB  NFLX   AAL   DAL   BAC  NVDA   WBD 
## 0.459 0.744 0.534 0.853 0.411 0.617 0.626 0.230 0.193 0.502 0.553 0.338 
## 
## Loadings:
##      Factor1 Factor2 Factor3
## MSFT  0.673   0.208   0.210 
## TSLA  0.462           0.190 
## AAPL  0.639   0.201   0.130 
## TWTR  0.148           0.348 
## AMZN  0.752           0.150 
## FB    0.568           0.237 
## NFLX  0.598           0.122 
## AAL   0.148   0.855   0.130 
## DAL           0.880   0.171 
## BAC   0.298   0.488   0.414 
## NVDA  0.634   0.201         
## WBD   0.101   0.220   0.776 
## 
##                Factor1 Factor2 Factor3
## SS loadings      2.870   1.934   1.135
## Proportion Var   0.239   0.161   0.095
## Cumulative Var   0.239   0.400   0.495
## 
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 40.08 on 33 degrees of freedom.
## The p-value is 0.185

6.Risk Management

Return <- return
Price <- asset
n = nrow(Return)
alpha = 0.05
# Non-paramatric

VaR_non = c()
ES_non = c()
for (i in 2:13) {
  q = as.numeric(quantile(unlist(Return[,i]), alpha)) #-0.06918179
  VaR_nonp = -100000 * q
  IEVaR = (Return[,i] < q)
  ES_nonp = -100000 * sum(Return[,i] * IEVaR) / sum(IEVaR)
  options(digits = 5)
  VaR_non = c(VaR_non, VaR_nonp)
  ES_non = c(ES_non, ES_nonp)
}
VaR_non

##  [1]  6918.2 15147.2 10749.9 22137.6  8881.9 10561.0 13072.7 15706.1 12891.9
## [10] 13041.2 11974.8 14903.9

ES_non

##  [1]  9183 20688 13639 25981 12875 15968 20587 24467 20512 16550 21517 19432

# Parametric - normal dist

# m = mean(Return$MSFT)
# sd = sd(Return$MSFT)
# qalpha = qnorm(alpha)
# VaR_par = -100000 * (m + sd * qalpha)
# ES = 100000 * (-m + sd * (dnorm(qalpha))/alpha)

VaR = c()
ES = c()
for (i in 2:13) {
  m = mean(unlist(Return[,i]))
  sd = sd(unlist(Return[,i]))
  qalpha = qnorm(alpha)
  VaR_par = -100000 * (m + sd * qalpha)
  ES_par = 100000 * (-m + sd * (dnorm(qalpha))/alpha)
  VaR = c(VaR, VaR_par)
  ES = c(ES, ES_par)
}
VaR

##  [1]  7239.2 23648.3  9963.7 23401.8 10831.8 11724.8 16743.7 18793.0 14803.4
## [10] 12075.0 14227.4 18533.6

ES

##  [1]  9712 30876 13200 29475 14243 15119 21650 23596 18777 15498 19149 23269

Bootstrap

library(bootstrap)

## 
## Attaching package: 'bootstrap'

## The following object is masked from 'package:pls':
## 
##     crossval

b=1000; n=400; S0=100000; var.boot=rep(0,b); es.boot=rep(0,b)

#non-parametric bootstrapping VaR
CI_nonp = c()
ES_nonp = c()
for (k in 2:13){
  for (i in 1:b){
    data.boot=sample(unlist(Return[,k]),n, replace=TRUE)
    var.boot[i]= -S0 * q
    #ES
    q = as.numeric(quantile(data.boot, 0.05))
    IEVaR = (data.boot < q)
    es.boot[i]  = -100000 * sum(data.boot * IEVaR) / sum(IEVaR)
    options(digits = 5)
    }
    CI_nonp = rbind(CI_nonp, c(quantile(var.boot, 0.025), quantile(var.boot, 0.975)))
    ES_nonp = rbind(ES_nonp, c(quantile(es.boot, 0.025), quantile(es.boot, 0.975)))
}

CI_nonp

##          2.5%   97.5%
##  [1,]  6610.1  7642.1
##  [2,] 14710.9 20336.6
##  [3,]  8108.5 11679.8
##  [4,] 19101.8 26024.1
##  [5,]  7861.3 10283.0
##  [6,]  9596.2 11187.7
##  [7,] 12068.5 19338.2
##  [8,] 14919.7 19849.6
##  [9,] 10980.4 15905.7
## [10,] 11456.9 13239.5
## [11,]  8707.3 18231.9
## [12,] 13540.7 18046.4

ES_nonp

##          2.5% 97.5%
##  [1,]  8062.5 10881
##  [2,] 18381.3 22243
##  [3,] 12208.0 15841
##  [4,] 24251.5 27276
##  [5,] 10885.7 16104
##  [6,] 12315.7 22335
##  [7,] 17424.7 24924
##  [8,] 20042.6 30785
##  [9,] 16429.8 27491
## [10,] 14336.8 20078
## [11,] 17571.3 24602
## [12,] 17529.6 21865

#parametric
CI_p = c()
ES_p = c()
for (k in 2:13){
  for (i in 1:b){
    data.boot=sample(unlist(Return[,k]),n, replace=TRUE)
    var.boot[i]= -S0*(mean(data.boot)+qnorm(0.05)*sd(data.boot))
    #ES
    s = sd(data.boot)*sqrt((1000-1)/1000)
    m = mean(data.boot)
    es.boot[i] = 100000 * (-m + s * (dnorm(qalpha))/alpha)
    }
    sd = sd(var.boot)*sqrt((1000-1)/1000)
    mean = mean(var.boot)
    CI_p = rbind(CI_p, c(mean + sd*qnorm(0.025), mean + sd*qnorm(0.975)))
    ES_p = rbind(ES_p, c(quantile(es.boot, 0.025), quantile(es.boot, 0.975)))
}
CI_p

##          [,1]    [,2]
##  [1,]  6279.0  8100.5
##  [2,] 21236.0 25752.1
##  [3,]  8682.5 11105.3
##  [4,] 21286.4 25179.3
##  [5,]  9567.7 11902.3
##  [6,]  9846.9 13421.5
##  [7,] 14874.3 18447.2
##  [8,] 16834.6 20502.2
##  [9,] 12629.6 16745.7
## [10,] 10566.4 13464.4
## [11,] 12295.1 16018.5
## [12,] 17044.9 19856.0

ES_p

##          2.5% 97.5%
##  [1,]  8574.4 10711
##  [2,] 27918.9 33628
##  [3,] 11648.3 14480
##  [4,] 27096.8 31459
##  [5,] 12741.4 15430
##  [6,] 12972.6 17073
##  [7,] 19407.2 23767
##  [8,] 21294.1 25589
##  [9,] 16312.3 21192
## [10,] 13873.4 17146
## [11,] 16796.7 21127
## [12,] 21350.5 24862

7.Copulas

library(copula)

#12 assets
est.df = data.frame()
i = 2
while (i <= 13) {
  est = as.numeric( fitdistr(unlist(Return[,i]),"t")$estimate )
  est[2] = est[2] * sqrt( est[3] / (est[3]-2) )
  est.row = c(est[1], est[2], est[3])
  est.df = rbind(est.df, est.row)
   i = i + 1
}
colnames(est.df) = c('m', 'sd', 'df')
# est.df: dataframe for estimations
est.df

##              m       sd      df
## 1   0.02427852 0.059806  6.1502
## 2   0.02107033 0.179692  3.8043
## 3   0.02794368 0.077174 85.9747
## 4   0.00336264 0.144514 26.7937
## 5   0.02330902 0.081512  8.4090
## 6   0.01693417 0.080778  5.0603
## 7   0.02101915 0.118101  5.4579
## 8   0.00098099 0.115981 10.4759
## 9   0.00847228 0.094574  4.9511
## 10  0.01587257 0.081778  7.2435
## 11  0.05033777 0.117408  9.9174
## 12 -0.01123768 0.119099  3.9940

data1 = cbind(rank(Return$MSFT)/(n+1), rank(Return$TSLA)/(n+1), rank(Return$AAPL)/(n+1), rank(Return$TWTR)/(n+1), rank(Return$AMZN)/(n+1), rank(Return$FB)/(n+1), rank(Return$NFLX)/(n+1), rank(Return$AAL)/(n+1), rank(Return$DAL)/(n+1), rank(Return$BAC)/(n+1), rank(Return$NVDA)/(n+1), rank(Return$WBD)/(n+1))

#fit t copula
cop_t_dim12 = tCopula(dim = 12, dispstr = "un", df=n-1)
ft1 = fitCopula(data = data1, copula = cop_t_dim12, method="ml")
ft1@estimate

##  [1]  0.85475  0.91235  0.80902  0.90781  0.88447  0.87286  0.85582  0.84566
##  [9]  0.89127  0.89719  0.87230  0.87739  0.82898  0.86659  0.86716  0.86022
## [17]  0.81634  0.80362  0.84942  0.85955  0.84101  0.83140  0.89427  0.89283
## [25]  0.85483  0.85002  0.84778  0.85846  0.89800  0.86533  0.84412  0.84790
## [33]  0.83700  0.79633  0.80647  0.82372  0.81253  0.82629  0.89189  0.91085
## [41]  0.83263  0.82139  0.86491  0.89658  0.84579  0.88087  0.82927  0.81514
## [49]  0.87028  0.85706  0.84999  0.81481  0.79458  0.84251  0.87361  0.83042
## [57]  0.94989  0.89428  0.84669  0.85586  0.88559  0.83657  0.85692  0.84774
## [65]  0.89854  0.82264 52.47837

ft1@loglik

## [1] 1445.5

AIC(ft1)

## [1] -2757.1

# fit gaussian copula
fnorm = fitCopula(data = data1, copula = normalCopula(dim=12), method="ml")
fnorm@estimate

## [1] 0.85056

fnorm@loglik

## [1] 1345.8

AIC(fnorm)

## [1] -2689.6

# frank
ffrank = fitCopula(copula = frankCopula(3, dim = 12), data = data1, method="ml")
ffrank@estimate

## [1] 15.113

ffrank@loglik

## [1] 1395.4

AIC(ffrank)

## [1] -2788.7

# clayton
fclayton = fitCopula(copula = claytonCopula(3, dim = 12), data = data1, method="ml")
fclayton@estimate

## [1] 1.0806

fclayton@loglik

## [1] 1078.3

AIC(fclayton)

## [1] -2154.7

# gumbel
fgumbel = fitCopula(copula = gumbelCopula(3, dim = 12), data = data1, method="ml")
fgumbel@estimate

## [1] 3.4771

fgumbel@loglik

## [1] 1301

AIC(fgumbel)

## [1] -2600.1

# joe
fjoe = fitCopula(copula = joeCopula(2, dim = 12), data = data1, method="ml")
fjoe@estimate

## [1] 13.161

fjoe@loglik

## [1] 1420

AIC(fjoe)

## [1] -2838.1